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Abstract 

The formulation and the implementation of 
boundary conditions within the context of the 
quadrature- free form of the discontinuous Galerkin 
method are presented for several types of bound- 
ary conditions for the Euler equations. An impor- 
tant feature of the discontinuous Galerkin method 
is that the interior point algorithm is well behaved 
in the neighborhood of the boundary and requires 
no modifications. This feature leads to a simple 
and accurate treatment for wall boundary conditions 
and simple inflow and outflow boundary conditions. 
Curved walls are accurately treated with only minor 
changes to the implementation described in earlier 
work. The “perfectly matched layer” approach to 
nonreflecting boundary conditions is easily applied 
to the discontinuous Galerkin. The compactness of 
the discontinuous Galerkin method makes it better 
suited for buffer-zone-type methods than high-order 
finite-difference methods. Results are presented for 
wall, characteristic inflow and outflow, and nonre- 
flecting boundary conditions. 

Introduction 

Much of the recent work in computational aeroa- 
coustics (CAA) has focused on improvements to tra- 
ditional finite-difference methods to increase the ac- 
curacy and to implement specialized boundary con- 
ditions. While this approach has promoted a rapid 
growth of the field, these methods place constraints 
on the mesh smoothness that make their application 
to highly complex geometries problematic. Further- 
more, the improved spatial operators are not appli- 
cable in the neighborhood of some critical flow phe- 
nomenon, such as shock waves, with out substantial 
modifications. The goal of this work is to develop 
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robust and efficient methods that pve accurate so- 
lutions independent of grid smoothness. 

The discontinuous Galerkin method is a highly 
compact formulation that provides a method of ob- 
taining the high accuracy required for CAA on non- 
smooth unstructured grids. The ability to use an 
unstructured grid greatly simplifies the largest ob- 
stacle in computing the flow around complex geome- 
tries: the generation of the grid. In reference 1, the 
discontinuous Galerkin method w-ts formulated in 
a quadrature-free form that reduced the computa- 
tional effort and storage requirements. In that work, 
the method was described in detail along with basic 
benchmark cases that demonstrate ihe accuracy and 
robustness of the method for the scalar advection 
equation and for the linear Euler equations. That 
work focused on the new implementation of the in- 
terior point scheme and addressed only periodic do- 
mains. 

In this article, the formulation ind implementa- 
tion of several types of boundary conditions for the 
linear Euler equations are described. Also discussed 
are features of the discontinuous Galerkin method 
that make the application of boundary conditions 
relatively straightforward and robust. These ben- 
eficial features are all attributable to the inherent 
compactness of the discontinuous Galerkin method. 

Most methods used for CAA today fall in 
the category of high-order finite difference meth- 
ods such, as the widely used di persion-relation- 
preserving (DRP) scheme? Effort; to develop spe- 
cialized boundary conditions for problems particu- 
lar to aeroacoustics have focused on finite-difference 
methods, but much of the work is also applicable to 
the discontinuous Galerkin method. In some cases, 
such as the work on wall boundary conditions by 
Tam and Dong, 3 special boundary conditions are 
needed to counter errors associaied with the ap- 
plication of finite-difference methods near a bound- 
ary: errors that do not occur in ffie discontinuous 
Galerkin method. 

The most problematic boundary in CAA is the 
boundary that is produced when an infinite or semi- 
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infinite domain is truncated to a finite computa- 
tional domain. In this case, precise flow conditions 
are not known at the boundary of the computational 
domain, and the boundary condition becomes more 
of a goal than a precise mathematical statement. In 
particular, the boundary condition seeks to make 
the flow field behave as if the computational domain 
were larger; waves are allowed to exit the computa- 
tional domain with no nonphysical side affects. In 
CAA, the boundary conditions appropriate for this 
type of boundary are referred to a s nonreflecting 
boundary conditions . 

Boundary conditions used for steady and unsteady 
aerodynamic calculations have relied primarily on 
characteristic formulations, such as the simple rela- 
tions proposed by Jameson et al. 4 to ensure that the 
correct information enters and leaves the domain; 
however, these methods become less accurate as the 
size of the computational domain is reduced. Efforts 
to improve on this have taken many forms, which 
range from efforts to analytically solve a simplified 
equation in the infinite domain outside the computa- 
tional domain 5, 6 to methods that solve specialized 
equations at the boundary or in a small region near 
the boundary. 2, 7 ’ 9l 10, 11 The methods work well 
when an acoustic wave exits the domain normal to 
the boundary; however, in other cases these meth- 
ods produce predictable reflections that depend on 
the angle of incidence in a manner that is fairly well 
understood in most cases. 

Two exceptions to this are the perfectly matched 
layer (PML) method of Hu 10 and the asymp- 

Q 

totic method of Hagstrom and Goodrich and 
Hagstrom. 11 Hagstrom’s approach is similar to that 
of Engquist and Majda 7 Giles, 9 and many others, 
except that the use of a Pade series approximation 
leads to a convergent sequence of equations in which 
the error associated with the wave incidence is re- 
duced as more terms are retained. In the PML ap- 
proach, a split and damped form of the governing 
equations is solved in a finite region near the bound- 
ary. Under certain constraints, no reflection of a 
wave of any incidence occurs at the interface be- 
tween the main computational domain and the re- 
gion where the PML method is applied. Within the 
PML region, waves are damped such that any reflec- 
tion of the wave off the outer boundary of the PML 
zone is insignificant. Because of the compact na- 
ture of the discontinuous Galerkin method, the PML 
method is more easily implemented for the discon- 
tinuous Galerkin method than for finite-difference 
methods. 

The first section of this article briefly describes the 
discontinuous Galerkin method and the quadrature- 


free form of the implementation; the reader is re- 
ferred to reference 1 for complete details. The second 
section describes issues related to boundary condi- 
tions and outlines the general approach to applying 
boundary conditions. The remaining sections deal 
with special issues of curved-wall and nonreflecting 
boundary conditions. Treatment of curved walls re- 
quires a minor modification to the basic formulation. 
Two types of nonreflecting boundary conditions are 
presented: a simple characteristic approach and the 
PML method. 


Discontinuous Galerkin Method 

The discontinuous Galerkin method is applicable 
to systems of first-order equations of the form 

^ + VF(t/) = 0 (1) 

defined on some domain Ct with a boundary 
where U = {u 0 , . . .} and F = {/o> /i , . . .}. The 

domain is partitioned into a set of nonoverlapping 
elements that cover the domain Q = G 

V t 

Within each element, the following set of equations 
is solved: 

/ bk^-Jidn - j Vb k -3~ l F i (V)J i dil 
n, 0 n t 

+ J b k J- 1 F R J i ds = 0 (2) 

an, 

for k = 1, 2, . . . , A, where k = 1,2,..., N] is a 
set of basis functions, 


N 


u 


Vi — ^ v tjbj , 

j = i 


J. = 


d(x 1 y y z) 


and Ji = |J,*|. Equation set (2) is obtained by pro- 
jecting equation (1) onto each member of the basis 
set and then integrating by parts to obtain the weak 
conservation form. In the present work, the basis set 
is the set of polynomials that are defined local to the 
element and are of degree < n. In two dimensions, 
for example, the basis set is {1,£, t? } £ 2 ,£?7, g 2 , . . .}, 
where (£, 77 ) are the local coordinates. The solu- 
tion U is approximated as an expansion in terms of 
the basis functions; thus, both V and F are discon- 
tinuous at the boundary between adjacent elements 
(hence, the name discontinuous Galerkin). The dis- 
continuity in V between adjacent elements is treated 

with an approximate Riemann flux, which is denoted 

F K; 

J, is the Jacobian of the transformation from 
the global coordinates (x,y, z) to the element coor- 
dinates (£,77,C) °f element i. Research has shown 12 
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that the upwind bias provided by the Riemann flux 
plays an important role in ensuring the stability of 
the discontinuous Galerkin method. 

In the usual implementation of the discontinu- 
ous Galerkin method/ the integrals are evaluated 
with quadrature formulas. This approach is prob- 
lematic for even moderately high-order implementa- 
tions in multidimensions and has limited most ef- 
forts to n . = 2 or 3. The difficulties arise in part 
because the number of quadrature points in multidi- 
mensional formulas of the required accuracy usually 
exceeds N (the number of terms in the expansion) by 
a considerable margin. In the quadrature- free form, 
the integral evaluations are reduced to a summa- 
tion over the coefficients of the solution expansion, 
which is an operation of order N . To implement the 
quadrature- free approach, the flux F must also be 
written as an expansion in terms of basis functions: 

N 

F(U)*G(V i )='£,9jW) b J 

j = 1 

(a similar expansion is made for the approximate 
Riemann flux F R .) This procedure is trivially ac- 
complished for linear equations, such as those of in- 
terest here. Several approaches to treating nonlinear 
equations are discussed and demonstrated in refer- 
ence 1. With the additional assumption that J, and 
J, are constant within each element, the integrals 
can be evaluated exactly, and the equation set can 
be rewritten in matrix form as 

m, 

+ B,,fc (j<Jr l 6&) • *i.fc = o (3) 

* = 1 

where m; is the number of sides around element t, 
Si k is the outward unit normal on side k ) V, = 
[v»,o, v*\ii • . .] T , Gj — 5*,i i • • •] i and G- ^ 

[g^ k o >9i*k i» • • -] T - The mass matrix M, and the vec- 
tor matrix A, are NxN matrices given by 

Mi = [m*,/], 

mkj = / bkbi dQ, 
n f 

for 1 < k,l < N . 

Derivation of the boundary integral terms is com- 
plicated only by the fact that the solutions on ei- 
ther side of the element boundary are represented 
in terms of different coordinate systems. This prob- 
lem is resolved by expressing the solution on both 


A, = [d k> t 


5 k ,i — / biVbk dQ 
a, 


sides of the element boundary in te ms of a common 
edge-based coordinate system (a simple coordinate 
transformation). This allows the boundary integral 
to be expressed in terms of an edge matrix Bp* times 
a vector that is composed of the coefficients of the 
approximate Riemann flux expresse d in terms of the 
edge-based coordinate system G^ ; (instead of the 
local element coordinate system < s in the case of 
Gi). 

In addition to the requirement that J t and J[ 
be constants within the element, most elements are 
constrained to shapes that map into one of a few 
fixed simple computational elemenls (such as a unit 
square or an equilateral triangle in two dimensions). 
With this last constraint, the matrices M,, A; and 
B,- * are the same for all elements of a given type, 
and the products M” 1 A and M“ l B/c can be pre- 
computed and stored at a considerable savings in 
terms of both computer storage and computational 
time. This constraint is only to facilitate an efficient 
implementation and can be relaxed at selected ele- 
ments if the need arises (e.g., to treat curved walls). 
A detailed derivation of the matrices M, A, and 
B k is given in reference 1. The special case of ele- 
ments with curved sides is described in a later sec- 
tion. Because equation (3) is of the same form for 
all elements, the element index i w 11 be dropped for 
clarity. 

Equation (3) is advanced in time by using 
the three-stage Runge-Kutta method of Shu and 
Osher. 16 Analysis of the stability of this approach 
can be found in reference 1 

General Features of Boundan v Conditions 

The first two terms of equation (3) depend only on 
the solution within the element, and communication 
between adjacent elements occurs only through the 
Riemann flux Gt . An important feature of the dis- 
continuous Galerkin method is thai the approximate 
Riemann flux is the only mechanism through which 
an element communicates with it surroundings, re- 
gardless of whether the element boundary is in the 
interior of the domain or coincides with the domain 
boundary. A notable consequence is that the usual 
interior algorithm is valid in elements adjacent to 
the boundary. In contrast, the interior point oper- 
ator of most high-order finite-diff-rence and finite- 
volume methods cannot be applied at points near 
the boundary without some modifications. These 
modifications usually result in reduced accuracy, and 
careful attention is required to prevent the introduc- 
tion of instabilities. 17 Thus, by us< of the discontin- 
uous Galerkin method, a major so iree of error com- 
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mon to many high-order finite-difference and finite- 
volume methods is completely avoided. 

Because each element communicates with its 
neighbors only through the approximate Riemann 
flux, most boundary conditions will be imposed via 
the approximate Riemann flux. In this respect, the 
imposition of boundary conditions for the discon- 
tinuous Galerkin method is quite similar to that of 
low-order finite-volume methods. This similarity is 
especially true for the quadrature-based discontinu- 
ous Galerkin method, in which the approximate Rie- 
mann flux is evaluated at discrete boundary points 
and then numerically integrated. In the quadrature- 
free form of the discontinuous Galerkin method, the 
approximate Riemann flux is a polynomial function 
of the edge coordinate and is never evaluated at spe- 
cific points. Thus, boundary conditions are applied 
to each component of the flux polynomial, rather 
than to the flux at specific points. 

Boundary conditions can be imposed either by 
providing the exterior side of the Riemann flux with 
a complete solution or by reformulating the bound- 
ary flux subject to the specified boundary conditions 
such that only the interior data is needed. However, 
either approach can be expressed exactly in terms of 
the other when the equations are linear. The first 
approach seems trivial to implement; however, this 
approach has the drawback that in most cases the 
complete solution is not known. Instead, the com- 
plete exterior solution must be reconstructed from 
the given boundary condition data and the interior 
solution. The work involved in a carefully derived re- 
construction procedure is usually equivalent to the 
work required to evaluate a completely reformulated 
flux, although the use (or misuse) of simple extrap- 
olation formulas is common. In this work, the ap- 
proximate Riemann flux on the boundary is replaced 
by a reformulated boundary flux. 

In the following discussion and examples, the lin- 
ear Euler equations in two dimensions will be used: 



- p-p - 


' 0 ' 

u = 

p 

, F = MU + 
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u 


iP 


V 
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where M = [M X) M y ], V = [u, v], and i and j are the 
Cartesian unit vectors [1,0] and [0,1], respectively. 
The components of U are normalized perturbation 
quantities from a free-stream condition about which 
the linearization has been performed. The compo- 


nent of the flux normal to the boundary is given by 


F(U) = F(U)-n = M n U + 


0 

V n 

aP 

PP 


( 5 ) 


where M n = M ■ n, V n = V ■ n, a = i n y P = j • n, 
and n = J“ jT Jds is the boundary-normal vector for 
an arbitrary edge. 


Wall Boundary Conditions 

Wall boundary conditions correspond to the case 
in which M n = 0 and V n is specified. Both 
symmetry-plane and hard-wall boundary conditions 
state that no flow passes through the boundary; 
thus, V n — M n — 0. The symmetry-plane bound- 
ary condition should be simply a special case of a 
general, hard-wall boundary condition in which the 
wall is planar; however, most finite-difference and 
finite- volume methods must treat the two differently 
in order to obtain accurate results. With the dis- 
continuous Galerkin method, the treatment of the 
two is identical because the discontinuous Galerkin 
method is valid without modification in the element 
next to the boundary. 

A transpiration wall condition is one in which fluid 
passes through the boundary at a specified rate. An 
example is a flow in which blowing or suction is ap- 
plied to a surface. Another example that is relevant 
to aerocoustic applications occurs when a flow is 
separated into incident and scatter components and 
each component is simulated individually. Occasion- 
ally, the form of the incident wave is known exactly, 
so that only the scattered wave needs to be simu- 
lated. With these assumptions, the flux through the 
boundary is given by 


F(U) wall = 


o - 

V n 

aP 


PP 


( 6 ) 


The flux is evaluated by using the pressure from the 
interior element and a specified function for V n . The 
function for V n must be expressed as a polynomial 
of the edge coordinates. This expression can be ob- 
tained by either a Taylor’s expansion or a projection 
procedure. Because the solution within each element 
is a known polynomial function, the interior solution 
at the edge is always available without the use of ex- 
trapolation formulas. 

Figure 1 illustrates a simple application of wall 
boundary conditions. An acoustic pulse is generated 
by a pressure disturbance in the initial condition of 
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an otherwise undisturbed Mach 0.5 flow'. The ini- 
tial pressure disturbance is a Gaussian distribution 
with a half-width of 0.05, centered at (-0.25, 0.25) in 
the domain 0 < x,y < 1. A hard- wall condition is 
specified on the y — 0 boundary, and through-flow 
boundary conditions (to be discussed later) are spec- 
ified on the other three sides of the domain. The re- 
sults shown are for a discontinuous Galerkin method 
with n = 4 (fifth order) and with the domain parti- 
tioned by an 18 x 18 triangulated grid. At t = 0.4, 
the incident pulse has reached the lower wall and 
has produced a reflection. In this case, the hard-wall 
boundary is equivalent to a symmetry plane. Figure 
1(b) shows similar results in which the computation 
included the mirror image of the original compu- 
tational domain. The maximum difference between 
the solutions is less that 0.1 percent and is attributed 
to the treatment of the flux at y = 0. In the first case 
in which y = 0 is a wall, the flux at y = 0 is given 
by F(U ) = [0,0,0,/?P] T In the case for which the 
y = 0 line is within the domain, the flux is evaluated 
by using the Lax-Friedrichs flux as the approximate 
Riemann solver* 

F(U U} Ui) = [ F(U U ) + F(Ut) - A (U u - I/,)] /2 

where subscripts u and / denote the upper and lower 
sides of the flux and A is greater than or equal to the 
magnitude of the largest eigenvalue of Assum- 
ing that the solution above and below y = 0 evolve 
symmetrically, U u is the same as Ui except for the 
v component, which is an odd function of y. Thus, 
the flux at y = 0 becomes 

T (U u ,Ui) = [0,0,0, (3P U + Xv u ] 

Because of the symmetry of the solution and the 
convergence properties of the discontinuous Galerkin 
method, v u goes to zero at the rate of Ax' 1 **' 1 ; thus, 
both formulas are accurate representations of the 
flux and exhibit the expected convergence proper- 
ties as the mesh is refined. Note, however, that a 
low-order error may be introduced if the solution 
is not symmetric, and if wall boundary conditions 
are implemented by retaining the approximate Rie- 
mann flux and evaluating the exterior solution with 
a reflection of the interior solution (as is commonly 
done on low-order finite-volume methods.) The spe- 
cific form of the error depends on the form of the 
approximate Riemann solver. 

Conditions at Curved Walls 

Walls that are smoothly curved can be mod- 
eled with at least second-order accuracy by straight 
line segments. To improve the accuracy requires a 


few simple modifications to the im plementation de- 
scribed previously. The first change is to compute 
distinct matrices M -1 A and M _ , B^ for each ele- 
ment and each side of that element that lies on a 
curved boundary. The only other :hange is simply 
to recognize that the edge normal vector ds is now 
a polynomial function instead of a constant vector; 
thus, aP and (3P in equation (6) are products of 
polynomials. Illustrated for triangles in figure 2, 
a general triangle with one curved side is mapped 
(with constant Jacobian) to a simple regular trian- 
gle in which the deviation of one si le from its usual 
straight line path is approximated by a polynomial 
i7(() wall* Because the Jacobian J i: constant within 
the element, it can be taken outride the integral; 
thus all integrations, matrix inversions, and matrix 
multiplications can be done in advance of the simu- 
lation as in the usual implementat on. The primary 
overhead associated with a curved < lement is the ad- 
ditional storage required to store a distinct copy of 
the matrices for each curved element. 

Figure 3 shows two solutions in which an acoustic 
pulse has passed over a cylinder o produce a re- 
flection. In the extreme case show n, the cylinder is 
modeled with only two elements. In figure 3(b) the 
curved sides are approximated by cubic polynomials. 
In this test case, the cylinder has a radius of 1/2, 
and the incident pulse is produce 1 by a Gaussian 
pressure disturbance in the initial solution at x — 3, 
y = 0. This case is similar to problem 2 of Cate- 
gory I of the recent workshop The Second Workshop 
on Benchmark Problems for CAA except that the 
Gaussian half-width of the initial disturbance has 
been increased to 1.6 (8 times la ger) so that the 
incident pulse is well resolved on extremely coarse 
grids and the error is dominated b\ the resolution of 
the cylinder. Figure 4 gives the convergence history 
of the solution as the resolution m increased. The 
average length scale of an element is defined as 

j area of dom ain 
V number of elements 

and the error is measured relative o reference solu- 
tion computed on a fine grid (At = 0.0498). The 
error is defined as the the L 2 norm of the difference 
in pressure at a large number of points uniformly dis- 
tributed in the region 0.63 < r < 2.0, 0 < 0 < 7t/ 2. 
The case with the cubic approximation for the wall 
maintains a fifth-order rate of convergence over the 
range of grids tested. The rate of convergence for 
the case with the linear approxim ition for the wall 
drops to less then third order as the mesh is refined. 
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Inflow and Outflow Boundary Conditions 

Inflow, outflow, and nonreflective boundary con- 
ditions are often treated as different entities; how- 
ever, for any system of equations such as the Euler 
equations at subsonic conditions, characteristic in- 
formation simultaneously enters and leaves the do- 
main through these boundaries. Typically, inflow 
and outflow boundary conditions have concentrated 
on ensuring that the correct information enters the 
domain; nonreflective boundary conditions have con- 
centrated on ensuring that waves that are leaving the 
domain can do so as if the boundary did not exist. 

The simplest form of an inflow and outflow bound- 
ary condition is obtained by splitting the flux into 
characteristic components and grouping the compo- 
nents according to whether their associated wave 
is entering or leaving the domain. The splitting is 
given by 

F = PfAjP'H/ 

= P [A + ] P interior + P [A ] P ^exterior 

where P and [A] are the eigenvector matrix and the 
diagonal matrix composed of the eigenvalues of 
respectively, and [A^] are diagonal matrices com- 
posed of just the positive or negative elements of 
[A]. The exterior solution is usually set to zero; how- 
ever, the solution could be set to any desired value 
to accommodate the case in which a specified wave 
is to enter the domain. This approach has been used 
in the results shown thus far and provides a crude, 
nonreflective boundary condition in that waves that 
are nearly aligned with the boundary will exit with 
little reflection. The method of Thompson 19 is an 
analogous procedure formulated for finite-difference 
methods. 

The reason for the reflection is that when an out- 
going wave that is not aligned with the boundary is 
decomposed into boundary-normal and boundary- 
tangent characteristic components, the inbound 
boundary-normal characteristic component is not 
exactly 0. Yet in almost all characteristic-based 
boundary conditions the inbound boundary-normal 
component set set either to 0 or to some specified 
value that has no relation to any outgoing wave that 
might exist. Most attempts to improve the nonre- 
flective boundary condition involve derivation of a 
means to reconstruct an inbound boundary-normal 
characteristic contribution associated with outgoing 
waves. The most accurate of these methods^’ 9 ’ 11 
involves the solution of an additional partial differ- 
ential equation along the boundary. Thus far, these 
boundary conditions have only been formally de- 
rived for smooth (if not planar) boundaries for which 


the mean flow is strictly inflow or outflow over the 
entire boundary. 

Another approach, the finite-wave model, pro- 
vides a simple (algebraic) method for improving the 
accuracy in some cases. This boundary condition 
was developed to deal with nonlinear effects of the 
Euler equations; however, the method also accounts 
for wave orientation in a way that is applicable to 
the linear case. The linear analog of the finite-wave 
mode is a simple modification to the standard char- 
acteristic approach and will be referred to as the 
modified characteristic method. The directionality 
inherent in the usual characteristic splitting arises 
because the boundary flux is the flux in the direc- 
tion of the boundary normal. The direction associ- 
ated with the flux cannot be altered; however, char- 
acteristic decomposition could certainly be based on 
another direction. In fact, because the boundary of 
the domain may not have any relation to the sound 
produced within the domain, other directions should 
be considered for the characteristic decomposition. 
If a single identifiable acoustic source is assumed, 
then the finite-wave model performs a characteristic 
decomposition along the assumed path of the wave. 
The decomposition is obtained from the characteris- 
tic variables associated with the flux in a prescribed 
direction: 

F W (U) = F(U) - w = M W U + 


where w is a unit vector in a prescribed direc- 
tion and other quantities are defined as in equa- 
tion (5) with n replaced by uj. The solution at 
the boundary associated with waves that are leav- 
ing the domain in the direction of w is given by 
Ut = P w [I + ] Pw 1 ^interior> wliere Pw are the eigen- 
vectors of and [/ + ] is a diagonal matrix with 
elements equal to 1 if the corresponding eigenvalue 
of is positive and equal to 0 otherwise. The flux 
through the boundary is given by evaluating equa- 
tion (5) with the solution f/&. 

The standard and modified characteristic meth- 
ods are compared in figure 5. The test case is the 
cylinder problem previously described with the non- 
reflecting boundary conditions imposed at r w 5.3. 
At time t = 10 most of the physical waves have 
exited the domain, and the remaining disturbances 
are caused by unwanted reflections. The modified 
characteristic boundary condition has reduced the 
reflection to less than half that of the standard char- 
acteristic boundary condition but the general form 
of the reflection is unchanged. 


0 

v w 

a u )P 
PxuP 
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The PML Method 

The PML method is a buffer-zone technique 
that solves a modified set of equations in a region 
that surrounds the primary computational domain. 
The modified equations are obtained by splitting 
the equations in boundary-normal and boundary- 
tangent directions and adding low-order damping to 
the boundary normal equations. For example, at 
boundaries aligned with either i or j, 


dlh d i}-FW) 

dt dx 


-a x U\ 


(7) 


and 


m d {> ■ ?w) 

dt + dy 


-(TyUl 


( 8 ) 


where U = U\ + Ui . The damping coefficients a x and 
(T y must be chosen such that the component of a that 
is tangent to the boundary does not vary along the 
boundary. This condition leads to the constraint on 
a in corner regions illustrated in figure 6 . Research 
has shown that , 10 for the ideal case of plane waves 
and straight boundaries that intersect at right angles 
(i.e., rectangular domains), no reflection of acoustic 
or convective waves will occur at the interface be- 
tween the primary computational domain and the 
PML zone, regardless of the angle at which waves 
strike the interface. The underlying theory places 
no constraint on the variation of cr in the direction 
normal to the boundary, but in applications to finite- 
difference methods a must vary smoothly. In numer- 
ical tests by Hu 10 the boundary-normal component 
of <r was increased quadratically a s a function of the 
distance from the interface. 

When the PML method is applied to the discon- 
tinuous Galerkin method, cr does not need to be var- 
ied smoothly. Furthermore, using a constant value 
of cr throughout a PML zone is advantageous. In 
the present application of the PML method to the 
discontinuous Galerkin method, the equations are 
solved in a different, but equivalent, form. In PML 
zones that border on the physical domain, the sum 
of the two split equations is solved in combination 
with the boundary-tangent equation. For example, 
on a boundary where x — Constant, cr y = 0 and the 
equations can be rewritten as 


and 


?^+VF(U) = -ct x (U-U 2 ) (9) 

ot 


DU, 9 (> ' Wl) 
1 r+ a, - 


( 10 ) 


The first equation is the standard nterior operator 
modified only by a zeroth-order dissipation term; 
thus, this equation is easily implen ented within the 
existing program structure. In a corner region, the 
sum of the split equations is solve 1 in conjunction 
with either equation (7) or ( 8 ): 

^7 + V • F(U) = -<r x U + (cr,. - <x y ){/ 2 (11) 


and 


dlh 

dt dy 


= - (TyU 2 


( 12 ) 


Note, however, that if <r x = cr y = t instant through- 
out the corner region, then the individual compo- 
nents U\ or U 2 do not contribute o equation ( 11 ); 
thus, only equation ( 11 ) needs to be solved. 

In figures 7 and 8 , solutions obtained with the 
PML method are compared with those obtained 
with characteristic boundary conditions. The test 
problem is a square domain (—50 *: x, y < 50) with 
hard-wall boundary conditions applied on the top, 
bottom, and left boundaries and eh her a PML zone 
or a characteristic boundary condition applied at the 
right boundary. The unsteady flow is initiated by 
a unit Gaussian pressure disturbance with a half- 
width of 3, positioned at x = 25, y = 0. The pri- 
mary domain is partitioned with an 18x18 triangu- 
lated grid; the PML case has two layers of elements 
in which the values of cr are constant: cr x — 0.2 and 
a y — 0. The solutions are compared with a base- 
line case in which the right boundary is extended to 
x & 161. Figure 7 shows the solutions at t = 40, 
which is just after the initial pulse has reached the 
boundary. The solution obtained w ith characteristic 
boundary conditions has weak refiections, and the 
solution obtain with the PML method agrees well 
with the baseline. Figure 8 shows the solutions at 
a much later time (f = 180) when eflections off the 
solid walls have produced a complex wave pattern. 
The solution obtained by using the PML method 
still agrees well with the baseline solution, while 
the standard characteristic method shows additional 
features that can only be attributed to nonphysical 
reflections off the right boundary. 

Figure 9 shows the effect of inci easing the thick- 
ness of the PML layer (x 6 - 50) and varying the 
values of <r r . The error metric is the maximum de- 
viation of pressure from that of the baseline solu- 
tion measured along the line x — 48 for t < 200. 
The solid line denotes the case in which cr x was var- 
ied quadratically, as described in reference 10 ; the 
dashed lines denote cases in which a x is fixed at one 
of several values. Note that the d ta at xt = 50 re- 
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suit from use of the standard characteristic bound- 
ary condition. 

Finally, the PML method is applied to the cylinder 
problem show earlier in the region r > 5.3. Figure 
10 shows the maximum pressure difference from the 
baseline solution measured at r = 5, t < 12. In this 
case, the PML method is implemented by assigning 
a normal direction to each element in the PML zone. 
Even though the boundary is curved and the normal 
direction varies in each element, the PML method 
shows a considerable improvement over the modified 
characteristic method (r^ & 5.3). 

Concluding Remarks 

The application of several types of boundary con- 
ditions for the discontinuous Galerkin method is pre- 
sented. Because of the compact form of the method, 
the discontinuous Galerkin method is applicable 
near boundaries without modification; this feature 
eliminates a major difficulty encountered by most 
high-order methods. As a consequence, boundary 
conditions such as symmetry-plane, curved-wall, and 
characteristic inflow outflow, are easy to implement 
and highly accurate. With modified characteristic 
boundary conditions that account for the direction 
of wave propagation, reflections are reduced to about 
half that of the standard characteristic method. The 
perfectly matched layer (PML) method is easily ap- 
plied to the discontinuous Galerkin method. The 
discontinuous Galerkin method allows the damping 
parameters to be abruptly turned on and then held 
constant within the PML zone. Reflections are re- 
duced by an order of magnitude below that of char- 
acteristic boundary conditions, even in cases with 
curved boundaries. 
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b. Computed mirror image of primary domain. 
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Figure 1. Reflection of cylindrical pressure pulse off 
a flat wall compared with direct o mputation of pri- 
mary domain plus mirror image. 
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Figure 6. PML zones with consistent, values of a. 



b. Characteristic boundary condition. 



c. PML zone. 


Figure 7. Comparison of pressure with different 
treatment of right boundary: t = 0. 
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Figure 9. Effect of thickness of PML layer and value 
of cr for flow shown in figures 7 and 8. 


Figure 8. Comparison of pressure with different 
treatment of right boundary: t — 180. 


c. PML zone. 


Figure 10. Effect of thickness of PML layer and value 
of cr for flow shown in figure 5. 


b. Characteristic boundary condition 



